f = @(x) tanh(2*x);
%f = @(x) x.^2;
sigma = 0.25;

N=200; xs = linspace(-4,4,N); ys = linspace(-4,4,N); x=[]; y=[]; px=[]; pyx=[]; for i=1:N; for j=1:N; x(i,j) = xs(i); y(i,j) = ys(j); px(i,j) = normpdf(x(i,j),0,1); pyx(i,j) = normpdf(y(i,j),f(x(i,j)),sigma); end; end;
p = px .* pyx;
px = (sum(p,2) * ((xs(end)-xs(1))/N)) * ones(1,N);
py = ones(N,1) * (sum(p,1) * ((ys(end)-ys(1))/N));
pxy = p ./ (py + 0.001);
pxdoy = px;

subplot(2,3,1);
contourf(x,y,p,7);
xlim([-2,2]);
ylim([-2,2]);
set(gca,'DataAspectRatio',[1 1 1]);
title('p(x,y)');
xlabel('x');
ylabel('y');
subplot(2,3,2);
contourf(x,y,pxy,7);
xlim([-2,2]);
ylim([-2,2]);
set(gca,'DataAspectRatio',[1 1 1]);
title('p(x | y)');
xlabel('x');
ylabel('y');
subplot(2,3,3);
contourf(x,y,pyx,7);
xlim([-2,2]);
ylim([-2,2]);
set(gca,'DataAspectRatio',[1 1 1]);
title('p(y | x)');
xlabel('x');
ylabel('y');
subplot(2,3,4);
M=10000;X=randn(M,1);Y=f(X)+sigma*randn(M,1); scatter(X,Y,'.'); 
xlim([-2,2]);
ylim([-2,2]);
set(gca,'DataAspectRatio',[1 1 1]);
title('data');
xlabel('x');
ylabel('y');
subplot(2,3,5);
contourf(x,y,px,7);
xlim([-2,2]);
ylim([-2,2]);
set(gca,'DataAspectRatio',[1 1 1]);
title('p(x | do(y))');
xlabel('x');
ylabel('y');
subplot(2,3,6);
contourf(x,y,pyx,7);
xlim([-2,2]);
ylim([-2,2]);
set(gca,'DataAspectRatio',[1 1 1]);
title('p(y | do(x))');
xlabel('x');
ylabel('y');

print('-depsc2','an.eps');
